Perturbation Trasmission and Graph visualization¶
In [ ]:
# Import necessary libraries
import networkx as nx # For creating and manipulating graphs
import numpy as np # For numerical operations
import matplotlib.pyplot as plt # For plotting
from matplotlib.patches import Patch # For creating custom legend patches
# Set the plotting style to 'ggplot'
plt.style.use('ggplot')
# Set the default figure size to [10, 10]
plt.rcParams['figure.figsize'] = [10, 10]
SIRS model simulation¶
In [ ]:
def dW(delta_t: float) -> float:
"""Sample a random number at each call."""
return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))
def Euler_Maruyama_method(parameters):
"""
Perform Euler-Maruyama method to simulate an epidemiological model.
Parameters:
- parameters (dict): Dictionary containing model parameters.
Expected keys:
't_in': Start time of simulation.
't_end': End time of simulation.
'N': Number of time steps.
'mu': Natural birth and death rate.
'b0': Transmission parameter.
'b1': Transmission parameter.
'phi': Transmission parameter.
'gamma': Recovery rate.
'ni': Disease-induced death rate.
'alpha': Scaling factor for stochastic term.
'S_in': Initial susceptible population.
'I_in': Initial infectious population.
'R_in': Initial recovered population.
Returns:
- TS (numpy.ndarray): Array of time steps.
- Ss (numpy.ndarray): Array of susceptible population over time.
- Is (numpy.ndarray): Array of infectious population over time.
- Rs (numpy.ndarray): Array of recovered population over time.
"""
t_in = parameters['t_in']
t_end = parameters['t_end']
N = parameters['N']
mu = parameters['mu']
b0 = parameters['b0']
b1 = parameters['b1']
phi = parameters['phi']
gamma = parameters['gamma']
ni = parameters['ni']
alpha = parameters['alpha']
S_in = parameters['S_in']
I_in = parameters['I_in']
R_in = parameters['R_in']
# Calculate time step size
dt = float((t_end - t_in) / N)
# Create array of time steps
TS = np.arange(t_in, t_end + dt, dt)
assert TS.size == N + 1
# Initialize arrays for S, I, R populations
Ss = np.zeros(TS.size)
Is = np.zeros(TS.size)
Rs = np.zeros(TS.size)
# Set initial values
Ss[0] = S_in
Is[0] = I_in
Rs[0] = R_in
# Perform Euler-Maruyama method to simulate SIR model
for i in range(1, TS.size):
t = t_in + (i - 1) * dt
S = Ss[i - 1]
I = Is[i - 1]
R = Rs[i - 1]
# Calculate effective transmission rate with stochastic term
b0_tilde = b0 + alpha * dW(dt)
beta = b0_tilde * (1 + b1 * np.cos(2 * np.pi * t + phi))
# Calculate stochastic term
dW_dt = dW(dt)
# Update S, I, R populations using differential equations
Ss[i] = S + ((mu - mu * S - beta * S * I + gamma * R) * dt - (alpha / b0_tilde) * beta * S * I * dW_dt)
Is[i] = I + ((beta * S * I - ni * I - mu * I) * dt - (alpha / b0_tilde) * beta * S * I * dW_dt)
Rs[i] = R + (ni * I - mu * R - gamma * R) * dt
return TS, Ss, Is, Rs
# Example parameters dictionary
parameters = {
't_in': 0,
't_end': 1,
'N': 5000,
'mu': 0.009,
'b0': 36.4,
'b1': 0.38,
'phi': 1.07,
'gamma': 1.8,
'ni': 36,
'alpha': 0.25,
'S_in': 0.998,
'I_in': 0.002,
'R_in': 0.0
}
# Run simulation using Euler-Maruyama method
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)
# Calculate basic reproduction number R0
def calculate_R0(parameters):
"""
Calculate the basic reproduction number R0 for an SIR-type epidemiological model.
Parameters:
- parameters (dict): Dictionary containing model parameters.
Expected keys:
'mu': Natural birth and death rate.
'gamma': Recovery rate.
'ni': Disease-induced death rate.
'b0': Transmission parameter.
'phi': Transmission parameter.
'b1': Transmission parameter.
Returns:
- R0 (float): Basic reproduction number.
"""
mu = parameters['mu']
gamma = parameters['gamma']
ni = parameters['ni']
beta = parameters['b0'] * parameters['phi'] * parameters['b1']
R0 = beta / (mu + ni + gamma)
return R0
# Calculate R0
R0 = calculate_R0(parameters)
# Print basic reproduction number (R0)
print(f"Basic Reproduction Number (R0): {round(R0, 1)}")
# Interpretation of R0
if round(R0) < 1:
print("The infection will likely die out over time.")
elif round(R0) == 1:
print("The infection will remain stable in the population but not cause an outbreak.")
else:
print("The infection will likely spread and cause an epidemic.")
# Additional statistics
peak_infection = max(Is)
final_size = Ss[-1] + Is[-1] + Rs[-1]
print(f"Peak Infection: {peak_infection}")
print(f"Final Size of the Epidemic: {final_size}")
Basic Reproduction Number (R0): 0.4 The infection will likely die out over time. Peak Infection: 0.0030332320384369007 Final Size of the Epidemic: 1.0009432020818723
SIRS PLOTS¶
In [ ]:
# Function to create a plot
def plot_population(TS, data, label, color):
plt.figure(figsize=(10, 6))
plt.plot(TS, data, label=label, color=color)
plt.xlabel('Time')
plt.ylabel('Proportion')
plt.title(f'{label} Population Over Time')
plt.legend()
plt.grid(True)
plt.show()
# Plot Susceptible over time
plot_population(TS, Ss, 'Susceptible', 'blue')
# Plot Infected over time
plot_population(TS, Is, 'Infected', 'red')
# Plot Recovered over time
plot_population(TS, Rs, 'Recovered', 'green')
# Plot all three (S, I, R) on the same graph
plt.figure(figsize=(10, 6))
plt.plot(TS, Ss, label='Susceptible', color='blue')
plt.plot(TS, Is, label='Infected', color='red')
plt.plot(TS, Rs, label='Recovered', color='green')
plt.xlabel('Time')
plt.ylabel('Proportion')
plt.title('SIR Model Over Time')
plt.legend()
plt.grid(True)
plt.show()
NetworkX¶
1. Small-World Networks (Watts-Strogatz Model)¶
Small-world networks have high clustering like regular networks but also have short average path lengths between nodes. The Watts-Strogatz model can be used to generate small-world networks.
In [ ]:
np.random.seed(42)
# Generate a Watts-Strogatz small-world network
n = 300
k = 4 # Each node is connected to k nearest neighbors
p = 0.15 # Probability of rewiring each edge
G = nx.watts_strogatz_graph(n, k, p)
# Ensure there is at least one initially infected node
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
while not any(state == 'I' for state in initial_states):
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
nx.set_node_attributes(G, dict(zip(G.nodes, initial_states)), 'state')
# Define colors for each state and create a custom legend
state_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
legend_elements = [Patch(facecolor=color, edgecolor=color, label=state) for state, color in state_colors.items()]
# Fixed layout for consistent node positions in plots
pos = nx.spring_layout(G)
# Function to update node states based on SIRS dynamics
def update_states(G, t, Is, Rs, Ss):
for node in G.nodes:
state = G.nodes[node]['state']
if state == 'S' and np.random.random() < Is[t]:
G.nodes[node]['state'] = 'I'
elif state == 'I' and np.random.random() < Rs[t]:
G.nodes[node]['state'] = 'R'
elif state == 'R' and np.random.random() < Ss[t]:
G.nodes[node]['state'] = 'S'
# Simulation and plotting
infection_alive = True # Flag to track if infection is still spreading
# Euler-Maruyama method simulation
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)
# Determine plot intervals
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)
# Simulation and plotting
for t_index in plot_intervals:
# Update node colors based on their states
node_color = [state_colors[G.nodes[node]['state']] for node in G.nodes]
# Draw the graph with updated node colors
plt.figure(figsize=(10, 6))
nx.draw(G, pos, with_labels=False, node_color=node_color, node_size=50)
plt.legend(handles=legend_elements, loc='best')
plt.title(f'Small-World Network (Watts-Strogatz Model) at t={TS[t_index]:.2f}')
plt.show()
# Determine the range for the next state update
next_t_index = plot_intervals[np.searchsorted(plot_intervals, t_index) + 1] if t_index != plot_intervals[-1] else len(TS)
# Update states of nodes in the specified range
for t in range(t_index, next_t_index):
update_states(G, t, Is, Rs, Ss)
# Check if all nodes are in susceptible state ('S')
if all(G.nodes[node]['state'] == 'S' for node in G.nodes):
infection_alive = False
# If infection has died out, print message and break out of loop
if not infection_alive:
print("The infection has died out.")
break
# Calculate effective reproduction number Rt over time
Rt = calculate_R0(parameters) * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")
The infection has died out. Effective Reproduction Number (Rt) at the beginning: 0.0013022154795595405
Various metrics to analyze the graph's properties¶
In [ ]:
# Degree distribution
degree_sequence = [d for n, d in G.degree()]
# Average degree
avg_degree = sum(degree_sequence) / len(degree_sequence)
print(f"Average degree: {avg_degree}")
# Clustering coefficient
avg_clustering = nx.average_clustering(G)
print(f"Average clustering coefficient: {avg_clustering}")
# Average shortest path length
avg_shortest_path_length = nx.average_shortest_path_length(G)
print(f"Average shortest path length: {avg_shortest_path_length}")
# Diameter (may need to handle exceptions if graph is not connected)
try:
diameter = nx.diameter(G)
print(f"Diameter: {diameter}")
except nx.NetworkXError:
print("Graph is not connected.")
# Degree centrality
degree_centrality = nx.degree_centrality(G)
# Example: Print degree centrality of the first few nodes
print("\nDegree Centrality\n")
for node in list(G.nodes())[:5]:
print(f"Node {node}: Degree Centrality = {degree_centrality[node]}")
# Betweenness centrality
betweenness_centrality = nx.betweenness_centrality(G)
# Example: Print betweenness centrality of the first few nodes
print("\nBetweenness Centrality\n")
for node in list(G.nodes())[:5]:
print(f"Node {node}: Betweenness Centrality = {betweenness_centrality[node]}")
# Assortativity
assortativity = nx.degree_assortativity_coefficient(G)
print(f"Degree Assortativity: {assortativity}")
# Global efficiency
global_efficiency = nx.global_efficiency(G)
print(f"Global Efficiency: {global_efficiency}")
# Local efficiency
local_efficiency = nx.local_efficiency(G)
print(f"Local Efficiency: {local_efficiency}")
# Degree histogram
degree_sequence = [d for n, d in G.degree()]
plt.hist(degree_sequence, bins='auto', alpha=0.75)
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution')
plt.show()
# Community detection using Louvain method
import community
partition = community.best_partition(G)
print(f"Number of communities detected: {len(set(partition.values()))}")
# Community detection using Louvain method
partition = community.best_partition(G)
print(f"Number of communities detected: {len(set(partition.values()))}")
# Modularity
modularity = community.modularity(partition, G)
print(f"Modularity: {modularity}")
# Visualization of communities (optional)
pos = nx.spring_layout(G)
plt.figure(figsize=(12, 12))
nx.draw_networkx_nodes(G, pos, node_size=20, cmap=plt.cm.RdYlBu, node_color=list(partition.values()))
nx.draw_networkx_edges(G, pos, alpha=0.4)
plt.title('Graph with Communities (Louvain method)')
plt.show()
Average degree: 4.0 Average clustering coefficient: 0.31297142857142857 Average shortest path length: 5.746570281124498 Diameter: 12 Degree Centrality Node 0: Degree Centrality = 0.020080321285140562 Node 1: Degree Centrality = 0.01606425702811245 Node 2: Degree Centrality = 0.020080321285140562 Node 3: Degree Centrality = 0.028112449799196783 Node 4: Degree Centrality = 0.020080321285140562 Betweenness Centrality Node 0: Betweenness Centrality = 0.02705176631312766 Node 1: Betweenness Centrality = 0.017062091061331035 Node 2: Betweenness Centrality = 0.028341754057983386 Node 3: Betweenness Centrality = 0.07175942440050946 Node 4: Betweenness Centrality = 0.01647276093760317 Degree Assortativity: 0.029907628017542705 Global Efficiency: 0.20521358507622173 Local Efficiency: 0.3944952380952381
Number of communities detected: 16 Number of communities detected: 15 Modularity: 0.745494
2. Random graph (Erdős-Rényi)¶
In [ ]:
np.random.seed(42)
# Create an Erdős-Rényi random graph
G = nx.erdos_renyi_graph(n=200, p=0.15) # n: number of nodes, p: probability of edge creation
# Assign initial states to nodes based on S_in, I_in, R_in
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
for node, state in zip(G.nodes, initial_states):
G.nodes[node]['state'] = state
# Simulate SIRS dynamics on the graph
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)
pos = nx.spring_layout(G) # Fixed layout for all iterations
# Determine the indices of the time points to plot
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)
# Define colors for each state
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
# Create custom legend
legend_elements = [
Patch(facecolor='blue', edgecolor='blue', label='Susceptible (S)'),
Patch(facecolor='red', edgecolor='red', label='Infected (I)'),
Patch(facecolor='green', edgecolor='green', label='Recovered (R)')
]
for t_index in plot_intervals:
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
node_color = [node_colors[G.nodes[node]['state']] for node in G.nodes]
# Draw the graph with updated node colors
nx.draw(G, pos, with_labels=True, node_color=node_color)
plt.legend(handles=legend_elements, loc='best')
plt.title(f'Random graph (Erdős-Rényi) at t={TS[t_index]:.2f}')
plt.show()
# Update states of nodes based on SIRS simulation results up to the next plot interval
for t in range(t_index, len(TS) if t_index == plot_intervals[-1] else plot_intervals[plot_intervals.tolist().index(t_index) + 1]):
for node in G.nodes:
if G.nodes[node]['state'] == 'S':
if np.random.random() < Is[t]: # Infection dynamics
G.nodes[node]['state'] = 'I'
elif G.nodes[node]['state'] == 'I':
if np.random.random() < Rs[t]: # Recovery dynamics
G.nodes[node]['state'] = 'R'
elif G.nodes[node]['state'] == 'R':
if np.random.random() < Ss[t]: # Loss of immunity dynamics
G.nodes[node]['state'] = 'S'
# Calculate effective reproduction number Rt over time
Rt = R0 * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")
Effective Reproduction Number (Rt) at the beginning: 0.005045444444444444
3. Scale-Free Networks (Barabási-Albert Model)¶
Scale-free networks are characterized by a small number of highly connected nodes (hubs) and many nodes with only a few connections. The Barabási-Albert model is a popular method to generate scale-free networks.
In [ ]:
np.random.seed(42)
# Generate a Barabási-Albert scale-free network
n = 200
m = 3 # Number of edges to attach from a new node to existing nodes
G = nx.barabasi_albert_graph(n, m)
# Assign initial states to nodes based on S_in, I_in, R_in
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
for node, state in zip(G.nodes, initial_states):
G.nodes[node]['state'] = state
# Simulate SIRS dynamics on the graph
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)
# Simulate SIRS dynamics on the graph
pos = nx.spring_layout(G) # Fixed layout for all iterations
# Determine the indices of the time points to plot
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)
# Define colors for each state
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
# Create custom legend
legend_elements = [
Patch(facecolor='blue', edgecolor='blue', label='Susceptible (S)'),
Patch(facecolor='red', edgecolor='red', label='Infected (I)'),
Patch(facecolor='green', edgecolor='green', label='Recovered (R)')
]
for t_index in plot_intervals:
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
node_color = [node_colors[G.nodes[node]['state']] for node in G.nodes]
# Draw the graph with updated node colors
nx.draw(G, pos, with_labels=True, node_color=node_color)
plt.legend(handles=legend_elements, loc='best')
plt.title(f'Scale-Free Networks (Barabási-Albert Model) at t={TS[t_index]:.2f}')
plt.show()
# Update states of nodes based on SIRS simulation results up to the next plot interval
for t in range(t_index, len(TS) if t_index == plot_intervals[-1] else plot_intervals[plot_intervals.tolist().index(t_index) + 1]):
for node in G.nodes:
if G.nodes[node]['state'] == 'S':
if np.random.random() < Is[t]: # Infection dynamics
G.nodes[node]['state'] = 'I'
elif G.nodes[node]['state'] == 'I':
if np.random.random() < Rs[t]: # Recovery dynamics
G.nodes[node]['state'] = 'R'
elif G.nodes[node]['state'] == 'R':
if np.random.random() < Ss[t]: # Loss of immunity dynamics
G.nodes[node]['state'] = 'S'
# Calculate effective reproduction number Rt over time
Rt = R0 * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")
Effective Reproduction Number (Rt) at the beginning: 0.005045444444444444